---
title: 'Replication: "Traditional Leaders and the 2014-2015 Ebola Epidemic"'
author: "Authors: Peter van der Windt & Maarten Voors"
date: "June, 2019"
output:
  html_document:
    number_sections: true
    theme: united
    toc: true
    toc_depth: 3
  pdf_document:
    number_sections: true
    toc: true
    toc_depth: 3
---


```{r, include=FALSE}
rm(list=ls(all=TRUE)) 
gc()
knitr::opts_chunk$set(echo = TRUE)
```


```{r, include=FALSE}
library(psych)
library(Hmisc)
library(dplyr)
library(knitr)
library(kableExtra)
library(stargazer)
require(MASS)
library(readstata13)

```



# Data Sources


```{r, include=FALSE}

# Set working directory and load data
setwd("C:/Users/Peter/Dropbox/RNE_institutions_health/03_analysis")
#setwd("/Users/maartenvoors/Dropbox/RNE_institutions_health/03_analysis/")
```

```{r, include=TRUE}

# Set working directory and load data

#setwd("")
d <- read.dta13("03_output/Dataset.dta")
```

"Dataset.dta" builds on multiple data sources, which are combined in "Prep_Data.do". All data sources and code to create "Dataset.dta" are available on Dataverse.

The following data sources were publicly available at the onset of the study (See online appendix for data source information for each variable):

* "pnas.1518587113.sd01.xlsx" and "pnas.1518587113.sd02.xlsx" from Fang et al (2016): https://www.pnas.org/content/suppl/2016/03/22/1518587113.DCSupplemental
* "chiefs_survey_clean_out.dta" from Acemoglu et al (2014). These were hand-coded from the annex of Acemoglu et al (2014), accessed via https://economics.mit.edu/files/9084 
* "nps-2007_cm_clean.dta" and "nps-2007_hh_clean.dta" from Bellows and Miguel (2009), which was accessed via https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/29272

In addition, we make use of the following data sources:

* "Core.dta" 149 chiefdoms, names of the 149 chiefdom to which all other data sources are merged (hand coded)
* "chiefdom_neighbors_matrix.dta", a 149x149 matrix indicating neighboring chiefdoms
* "SL_admin_chiefdom.shp", a shapefile used to create Figure B1
* "ruggedness_by_chiefdom_SL_30m_1km.xlsx", a measure for a chiefdom's ruggedness (created using ArcGIS)

# Prepare variables

Prepare the dependent variables:

```{r, include=TRUE}

# In logs
d$ebola_confirmed_ln <- log(d$ebola_confirmed+ sqrt(d$ebola_confirmed^2 + 1 ) )

## Dependent variables used in robustness tests:

# Days between symptom and test
d$time= ((d$conf_TimeOnsetTest * d$ebola_confirmed) + (d$sus_TimeOnsetTest * d$ebola_suspected)) / (d$ebola_confirmed + d$ebola_suspected)

# Total cases
d$ebola_total <- d$ebola_confirmed + d$ebola_suspected
d$ebola_total_ln <- log(d$ebola_total + sqrt(d$ebola_total^2 + 1))

# Not taking into account the first 3 months
d$ebola_confirmedSUB3_ln <- log(d$ebola_confirmedSUB3+ sqrt(d$ebola_confirmedSUB3^2 + 1 ) )


```


Next, we prepare the independent variable:

```{r, include=TRUE}

# Paramount Chief power (inverse number families)
table(d$fam_num)
d$authority = 1/d$fam_num

```

Create a matrix with total Ebola cases in a chiefdom's neighboring chiefdoms:

```{r, include=TRUE}

d <- d[order(d$ID),] 
d$cases_chiefdomX <- d$ebola_confirmed_ln
# Make the index from long to wide
d$x <- seq(dim(d)[1])
d<-reshape(d, v.names="cases_chiefdomX", idvar="x", timevar="ID", direction="wide")

for (i in 1:149){
  d[, paste("cases_chiefdom", i, sep="")]<-NA
}

for (i in 1:149){
   d[, which(names(d)=='cases_chiefdom1')+i-1] <- d[i, which(names(d)=='cases_chiefdomX.1')+i-1]
}

d<-d[, !grepl("cases_chiefdomX.", colnames(d))] 

## Multiply neighbors with Ebola incidence
for (i in 1:149){
  d[, paste("Z", i, sep="")]<-NA
}

for (j in 1:149){
  for (i in 1:149){
    temp_neighbor<-d[j, which(names(d)=='neighbor1')+i-1]
    temp_cases_chiefdom<-d[j, which(names(d)=='cases_chiefdom1')+i-1]
    d[j, which(names(d)=='Z1')+i-1]<- temp_neighbor * temp_cases_chiefdom
  }
}

# Wy_total: total cases in neighboring chiefdoms (W is neighbour matrix, y is cases vector)
d$Wy_total <- rowSums(d[,(which(names(d)=='Z1'): which(names(d)=='Z149')) ], na.rm=TRUE)

# Wy_average: average cases in neighboring chiefdoms
d$t_neig <- rowSums(d[,(which(names(d)=='neighbor1'): which(names(d)=='neighbor149')) ], na.rm=TRUE)  
d$Wy_average <- d$Wy_total/d$t_neig

# drop t_neig Z* neighbor* cases_chiefdom*
d<-d[, !grepl("neighbor", colnames(d))]   
d<-d[, !grepl("cases_chiefdom", colnames(d))]    
d<-d[, !grepl("Z", colnames(d))]    
d$t_neig<-NULL      

```

Standardize all variables:

```{r, include=TRUE}

d$Sebola_confirmed <- scale(d$ebola_confirmed)
d$Sebola_suspected <- scale(d$ebola_suspected)
d$Sebola_confirmed_ln <- scale(d$ebola_confirmed_ln)
d$Sebola_total_ln<- scale(d$ebola_total_ln)
d$Sebola_confirmedSUB3_ln <-scale(d$ebola_confirmedSUB3_ln)
d$Sfam_num <-scale(d$fam_num)
d$Sauthority <- scale(d$authority)
d$Sn_chiefs <- scale(d$n_chiefs)
d$Samalgamation <- scale(d$amalgamation)
d$Smining_1930 <- scale(d$mining_1930)
d$Scoast_dist <- scale(d$coast_dist)
d$Sriver_dist <- scale(d$river_dist)
d$Smitchell_km <- scale(d$mitchell_km)
d$Srail_dist <- scale(d$rail_dist)
d$Sdist_min <- scale(d$dist_min)
d$SWy_total <- scale(d$Wy_total)
d$Stime <- scale(d$time)
d$Sk29 <- scale(d$k29)
d$Sk31 <- scale(d$k31)
d$Sk39 <- scale(d$k39)
d$Sk40 <- scale(d$k40)
d$Sh1a <- scale(d$h1a)
d$Sh1b <- scale(d$h1b)
d$Sh1c <- scale(d$h1c)       
d$Spop <- scale(d$pop)

```


# Table 1: Descriptive Statistics

```{r, include=TRUE}

summ <- function(x) c(length(x), mean(x), sd(x), min(x), max(x))

table_summ1 <- sapply(list(d$ebola_confirmed, d$fam_num, d$authority, d$n_chiefs, d$amalgamation, d$mining_1930, d$coast_dist, d$river_dist, d$mitchell_km, d$rail_dist, d$dist_min, d$pop), summ) %>% t()

rownames(table_summ1) <- c("Ebola confirmed cases", 
                          "Number of ruling families",
                          "Paramount Chief power", 
                          "Chiefs recalled",
                          "Amalgamation",
                          "Mining permission", 
                          "Distance to coast",
                          "Distance to nearest river", 
                          "Distance to 1895 trade routes", 
                          "Distance to 1907 railroad",
                          "Minimum distance to major city",
                          "Population size in 2004")

kable(table_summ1, caption = "Table 1: Descriptive statistics", col.names =  c( "N", "Mean", "SD"," Min" ,"Max") , digits = 2 ) %>% kable_styling(bootstrap_options = "striped", full_width = T) %>% footnote(general = "Detailed variable descriptions are in online Appendix A.", general_title = "Note: ", footnote_as_chunk = T)

```

# Table 2: Paramount Chief Power and Ebola Cases

```{r, results ="asis", warning=FALSE, include=TRUE}

#1: basic
t2.model.1 <- lm(Sebola_confirmed_ln ~ Sauthority, data=d)

#2: 1 + ARR controls + district FE
t2.model.2 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

#3: 2 but using negative binomial
t2.model.3 <- glm.nb(ebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

#4: 2 + controlling for population
t2.model.4 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + Spop + factor(DISTRICT), data=d)

#5: 4 + Spatial matrix
t2.model.5 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + Spop + factor(DISTRICT) + SWy_total, data=d)

#6: 5 + >0 cases
t2.model.6 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + Spop + factor(DISTRICT) + SWy_total, data=subset(d, ebola_confirmed>0))

stargazer(t2.model.1, t2.model.2, t2.model.3, t2.model.4, t2.model.5, t2.model.6, 
          type = "html", 
          out.header=TRUE, 
          omit = c("DISTRICT", "Constant"), 
          omit.stat = c("adj.rsq", "theta", "f", "aic", "ser"), 
          title = "Table 2. Paramount Chief Power and Ebola Cases", 
          covariate.labels = c("Paramount Chief power", "Chiefs recalled", "Amalgamation", "Mining", "Coast", "Nearest river", "Trade route", "Railroad", "Major city", "Population size", "Ebola neighbors"), 
          column.labels =c("Basic", "FEs+Controls", "Neg. Bin.", "Population", "Spatial", "Exposure"), 
          dep.var.labels.include = FALSE, 
          dep.var.caption  = "Model Specification",
          model.names = FALSE, 
          add.lines = list(c("District FEs", "N", "Y", "Y", "Y", "Y", "Y"), 
                           c("Regression", "OLS", "OLS", "Neg. Bin.", "OLS", "OLS", "OLS")), 
          notes = c("Regressions at chiefdom level. All variables are standardized. Standard errors in parentheses. * p<0.10 ; ** p<0.05 ; *** p<0.01 (two-tailed)"),
          notes.align = "l", 
          notes.append = FALSE
          )



```


Chiefdoms with a one standard deviation (SD) "stronger" Paramount Chief are associated with a decrease of 0.18 SD log Ebola cases (p<0.05). In substantive terms, this result implies that adding one additional ruling family (i.e. less powerful Paramount Chief) leads to an increase of 22.6% Ebola cases in that Chiefdom. 

```{r, echo=TRUE, include=TRUE}

# Model 2 but, for ease of interpretation, dependent and independent variables are not standarized and independent variable is simply number of families

model.inter <- lm(ebola_confirmed_ln ~ fam_num + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

# summary(model.inter)

# effect of one fewer family
-1 * (exp(model.inter$coefficients[2]) - 1) * 100

```


# Table C1: Paramount Chief Power and Ebola Cases, Terrain Control

```{r, results ="asis", warning=FALSE, include=TRUE}

#table(d$ruggmean1km)

#1: basic
tx.model.1 <- lm(Sebola_confirmed_ln ~ Sauthority + ruggmean1km, data=d)

#2: 1 + ARR controls + district FE
tx.model.2 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT)  + ruggmean1km, data=d)

#3: 2 but using counts using a negative binomial
tx.model.3 <- glm.nb(ebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT) + ruggmean1km, data=d)

#4: 2 + controlling for population
tx.model.4 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + Spop + factor(DISTRICT) + ruggmean1km, data=d)

#5: 4 + Spatial matrix
tx.model.5 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + Spop  +factor(DISTRICT) + SWy_total + ruggmean1km, data=d)

#6: 5 + >0 cases
tx.model.6 <- lm(Sebola_confirmed_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + Spop  +factor(DISTRICT) + SWy_total + ruggmean1km, data=subset(d, ebola_confirmed>0))


stargazer(tx.model.1, tx.model.2, tx.model.3, tx.model.4, tx.model.5, tx.model.6, 
          type = "html", 
          out.header=TRUE, 
          omit = c("DISTRICT", "Constant"), 
          omit.stat = c("adj.rsq", "theta", "f", "aic", "ser"), 
          title = "Table C1: Paramount Chief Power and Ebola Cases, Terrain Control", 
          covariate.labels = c("Chief strength", "Chiefs recalled", "Amalgamation", "Mining", "Coast", "Nearest river", "Trade route", "Railroad", "Major city", "Population size", "Ebola neighbors", "Ruggedness"), 
          column.labels =c("Basic", "FEs+Controls", "Neg. Bin.", "Population", "Spatial", "Exposure"), 
          dep.var.labels.include = FALSE, 
          dep.var.caption  = "Model Specification",
          model.names = FALSE, 
          add.lines = list(c("District FEs", "N", "Y", "Y", "Y", "Y", "Y"), 
                           c("Regression", "OLS", "OLS", "Neg. Bin.", "OLS", "OLS", "OLS")), 
          notes = c("Regressions at chiefdom level. All variables are standardized. Standard errors in parentheses. * p<0.10 ; ** p<0.05 ; *** p<0.01 (two-tailed)"),
          notes.align = "l", 
          notes.append = FALSE
          )


```


Also, when we run a regression, akin to the "balance test" in Table 3 in Acemoglu, Reed and Robinson (2013), we find that ruggedness and chiefly power are unrelated, with a p-value equal to:

```{r, results ="asis", warning=FALSE, include=TRUE}

model1 <- lm(ruggmean1km ~ Sauthority + Sdist_min + Sn_chiefs + Samalgamation + factor(DISTRICT), data=d)

summary(model1)$coefficients[2,4] 
```

The full model is as follows:
```{r, results ="asis", warning=FALSE, include=TRUE}

stargazer(model1, 
          type = "html", 
          out.header=TRUE, 
          omit = c("DISTRICT", "Constant"),
          omit.stat = c("adj.rsq", "theta", "f", "aic", "ser"), 
          title = "Chiefly Power and Ruggedness", 
          covariate.labels = c("Chief strength", "Chiefs recalled", "Amalgamation"), 
          dep.var.labels.include = FALSE, 
          model.names = FALSE, 
          add.lines = list(c("District FEs", "Y"), 
                           c("Regression", "OLS")), 
          notes = c("Regressions at chiefdom level. All variables are standardized. Standard errors in parentheses. * p<0.10 ; ** p<0.05 ; *** p<0.01 (two-tailed)"),
          notes.align = "l", 
          notes.append = FALSE
          )


```


# Table D1: Paramount Chief Power and Per Capita Ebola Cases

```{r Table A4, results ="asis", warning=FALSE, include=TRUE}

d$ebol_ln <- log(d$ebola_confirmed + sqrt(d$ebola_confirmed^2 + 1))
d$pop_ln <- log(d$pop + sqrt(d$pop^2 + 1))
d$ebol_pop <- d$ebol_ln - d$pop_ln
d$Sebol_pop <- scale(d$ebol_pop) 

#1: basic
ta2.model.1 <- lm(Sebol_pop ~ Sauthority, data=d)

#2: 1 + ARR controls + district FE
ta2.model.2 <- lm(Sebol_pop ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

#3: 2 + Spatial matrix
ta2.model.3 <- lm(Sebol_pop ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT) + SWy_total, data=d)

#4: 3 + >0 cases
ta2.model.4 <- lm(Sebol_pop ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT) + SWy_total, data=subset(d, ebola_confirmed>0))

stargazer(ta2.model.1, ta2.model.2, ta2.model.3, ta2.model.4, 
          type = "html", 
          out.header=TRUE, 
          omit = c("DISTRICT", "Constant"),
          omit.stat = c("adj.rsq", "theta", "f", "aic", "ser"), 
          title = "Table D1: Paramount Chief Power and Per Capita Ebola Cases", 
          covariate.labels = c("Chief strength", "Chiefs recalled", "Amalgamation", "Mining", "Coast", "Nearest river", "Trade route", "Railroad", "Major city", "Ebola neighbors"), 
          column.labels =c("Basic", "FEs+Controls", "Spatial", "Exposure"), 
          dep.var.labels.include = FALSE, 
          dep.var.caption  = "Model Specification",
          model.names = FALSE, 
          add.lines = list(c("District FEs", "N", "Y", "Y", "Y"), 
                           c("Regression", "OLS", "OLS", "OLS", "OLS")), 
          notes = c("Regressions at chiefdom level. All variables are standardized. Standard errors in parentheses. * p<0.10 ; ** p<0.05 ; *** p<0.01 (two-tailed)"),
          notes.align = "l", 
          notes.append = FALSE
          )


```


# Table E1: Additional Descriptive Statistics

```{r, results ="asis", warning=FALSE, include=TRUE}

summ2 <- function(x) c(length(na.omit(x)), mean(x, na.rm=TRUE), sd(x, na.rm=TRUE), min(x, na.rm=TRUE), max(x, na.rm=TRUE)) 

table_summ2 <- sapply(list(d$k29, d$k39, d$k31, d$time, d$k40, d$h1a, d$h1b, d$h1c, d$ebola_suspected, d$ebola_confirmedSUB3), summ2) %>% t()

rownames(table_summ2) <- c("Presence facility", 
                           "Presence doctor", 
                           "Presence nurse", 
                           "Time between symptom and test", 
                           "Cost of use", 
                           "Use government health facility",
                           "Use traditional health facility",
                           "Use pharmacist", 
                           "Ebola suspected cases", 
                           "Ebola confirmed cases (>3 months)")
                           
kable(table_summ2 , caption = "Table E1: Additional Descriptive Statistics", col.names =  c( "N", "Mean", "SD"," Min" ,"Max") , digits = 2 ) %>% kable_styling(bootstrap_options = "striped", full_width = T) %>% footnote(general = "Detailed variable descriptions are in the online appendix.", general_title = "Note: ", footnote_as_chunk = T)


```



# Table E2: Paramount Chief Power and Reporting of Ebola Cases

```{r, results ="asis", warning=FALSE, include=TRUE}

# Time between onset and test: combine confirmed and suspected	
ta4.model.1 <-lm(Stime ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

# Presence facility, etc.
ta4.model.2 <- lm(Sk29 ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d) 
ta4.model.3 <- lm(Sk39 ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d) 
ta4.model.4 <- lm(Sk31 ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d) 

# Cost to use, etc.
ta4.model.5 <- lm(Sk40 ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)
ta4.model.6 <- lm(Sh1a ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)
ta4.model.7 <- lm(Sh1b ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)
ta4.model.8 <- lm(Sh1c ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

# Total cases	
ta4.model.9 <- lm(Sebola_total_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

# Confirmed cases taking out the first few months	
ta4.model.10 <- lm(Sebola_confirmedSUB3_ln ~ Sauthority + Sn_chiefs + Samalgamation + Smining_1930 + Scoast_dist + Sriver_dist + Smitchell_km + Srail_dist + Sdist_min + factor(DISTRICT), data=d)

stargazer(ta4.model.1, ta4.model.2, ta4.model.3, ta4.model.4, ta4.model.5, ta4.model.6, ta4.model.7, ta4.model.8, ta4.model.9, ta4.model.10, 
          type = "html", 
          out.header=TRUE, 
          omit = c("DISTRICT", "Constant"),
          omit.stat = c("adj.rsq", "theta", "f", "aic", "ser"), 
          title = "Table E2: Paramount Chief Power and Reporting of Ebola Cases", 
          covariate.labels = c("Chief strength", "Chiefs recalled", "Amalgamation", "Mining", "Coast", "Nearest river", "Trade route", "Railroad", "Major city"), 
          column.labels =c("Process time", "Presence facility", "Presence doctor", "Presence nurse", "Cost to use", "Use gov.", "Use trade", "Use pharmacy", "Total", "Month >3"), 
          dep.var.labels.include = FALSE, 
          model.names = FALSE, 
          add.lines = list(c("District FEs", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y"), 
                           c("Regression", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS")), 
          notes = c("Regressions at chiefdom level. Regressions include the number of chiefs recalled and amalgamation of chiefdoms plus the six geographic controls. Standard errors in parentheses. * p<0.10 ; ** p<0.05 ; *** p<0.01 (two-tailed)"),
          notes.align = "l", 
          notes.append = FALSE
          )


```



End of replication file. 
